## ----setseedch3, echo = FALSE--------------------------------------------
set.seed(1982)


knitr::opts_chunk$set(
                      echo = FALSE,
                      warning = FALSE,
                      message = FALSE
)


## ----loadlibsch3---------------------------------------------------------
library(tidyverse)
library(reshape2)
library(ascii)
library(zoo)
library(lme4)
library(memisc)
library(pander)
library(Amelia)
library(arm)
library(broom)
library(gridExtra)
require(rstanarm)
library(stringr)
source("common_funcs.R")


## ----datain--------------------------------------------------------------
pta <- read.csv("./pta_data.csv")
pta$Date <- as.Date(pta$Date, "%d/%m/%Y")

### Make sure only to get decisions with a UKSC caseid
pta <- subset(pta, grepl("UKSC", pta$CASEID))

reports <- read.csv("./reports_data.csv")
reports <- unique(reports)
pta <- merge(pta, reports,
             by.x = c("Neutral"),
             by.y = c("APPORDNEUTRAL"),
             all.x = TRUE,
             all.y = FALSE,
             sort = FALSE)



## ----crossappeals--------------------------------------------------------
### Need to remove some cases w/ the same CASEID
### these are typically cross-appeals
pta <- pta %>% 
	group_by(CASEID) %>% 
    filter(Date == min(Date))
### Now, what about multiple applications?
pta$AppType <- factor(pta$AppType,
                      levels = LETTERS[1:10],
                      ordered = TRUE)
pta <- pta %>%
    group_by(CASEID) %>%
    filter(AppType == max(AppType))
pta$AppType <- as.character(pta$AppType)

### If we still need to remove observations
pta <- pta %>%
    group_by(CASEID) %>%
    filter(row_number() == 1)



## ----importance----------------------------------------------------------
pta$Reported <- gsub("Where Reported: ","",pta$Reported)
unique.reports <- lapply(strsplit(pta$Reported, ";"), function(x) gsub("[^A-Za-z ()]","",x))
unique.reports <- lapply(unique.reports, trim)
unique.reports <- lapply(unique.reports, function(x)
	gsub("Times [A-Z][a-z]+","Times",x))
unique.reports <- lapply(unique.reports, function(x)
	gsub("\\(\\) ?","",x))

kill.list <- c("Official Transcript", "NIQB", "NICA", "HCJAC",
	"EWHC (Ch)", "EWHC (QB)", "EWCA Civ", "EWCA Crim",
	"EWHC (Admin)", "CSIH") 
unique.reports <- lapply(unique.reports, function(x) 
	x[!is.element(x, kill.list)])

generalist <- c("All ER", "WLR", "SJLB", "Times", "SLT", "SC", "GWD", "NI")
generalist.count <- lapply(unique.reports, function(x) 
	x[is.element(x, generalist)])
specialist.count <- lapply(unique.reports, function(x) 
	x[!is.element(x, generalist)])

pta$Importance <- unlist(lapply(generalist.count, length))



## ----workload------------------------------------------------------------
load("./ch2_data.RData")
old <- nrow(pta)
pta <- merge(pta, combWorkload,
	by.x = "Date",
	by.y = "Date",
	all.x = T,
	all.y = F,
	sort = FALSE)
new <- nrow(pta)
stopifnot(old == new)
sheet1$PTAd <- sheet1$CASEID %in% unique(pta$CASEID)



## ----ptapercents---------------------------------------------------------
later.years <- which(sheet1$HANDDOWN > as.Date("2011-01-01"))
overall <- round(mean(sheet1$PTAd[later.years]) * 100)
ewni <- round(mean(sheet1$PTAd[later.years][which(sheet1$Area[later.years]!="Scots")]) * 100)
scots <- round(mean(sheet1$PTAd[later.years][which(sheet1$Area[later.years]=="Scots")]) * 100)


## ---- fig.cap = "Success rates of permission to appeal applications grouped by court releases. Solid line indicates a six-month rolling average. See text for further details. ", fig.scap = "Success rates of permission to appeal applications. ", fig.width = 6, fig.height = 4, dpi = 300, out.width = "100%"----
pta$success <- as.numeric(pta$Outcome %in% c("Granted","Granted in part"))

### Month by month plot of success rates
pta$yymm <- as.Date(paste0("20", format(pta$Date, "%y-%m"), "-01"))

plot.df <- pta %>%
    group_by(File) %>%
    summarize(nCases = n(),
              yymm = median(yymm, na.rm = TRUE),
              rate = mean(success, na.rm = TRUE)) %>%
    arrange(yymm) %>%
    mutate(rolling = rollmean(rate, 6, na.pad = TRUE))

ggplot(plot.df, aes(x = yymm, y = rate)) +
    scale_x_date("Date of PTA meeting") +
    scale_y_continuous("Proportion of cases granted permission",
                       labels = scales::percent) + 
    geom_point(colour = "grey50") +
    geom_line(aes(y = rolling)) +
    theme_uksc()



## ----ratesbyarea, fig = TRUE, fig.cap = "Success rates by area of law", fig.width = 7, fig.height = 4, dpi = 300, out.width = "100%"----
plot.df <- pta %>%
    group_by(Area) %>%
    summarize(nCases = n(),
              label = paste0(unique(Area), "\n(", nCases, " cases)"),
              rate = mean(success, na.rm = TRUE)) %>%
    filter(!is.na(Area))

plot.df$label <- factor(plot.df$label,
                       levels = plot.df$label,
                       ordered = TRUE)
plot.df$label <- reorder(plot.df$label, plot.df$rate)

ggplot(plot.df, aes(x = label, y = rate)) +
    scale_x_discrete("Area of law") +
    scale_y_continuous("Proportion of cases granted permission",
                       labels = scales::percent) + 
    geom_bar(stat = "identity", colour = "grey50") +
    theme_uksc()



## ----ratesbyimportance, fig = TRUE, fig.cap = "Success rates in cases according to their legal importance. ", fig.width = 9, figh.height = 4----
plot.df <- pta
plot.df$Importance[which(plot.df$Importance > 4)] <- 4
plot.df <- plot.df %>%
    group_by(Importance) %>%
    summarize(nCases = n(),
              label = paste0(unique(Importance), "\n(", nCases, " cases)"),
              rate = mean(success, na.rm = TRUE)) %>%
    filter(!is.na(Importance))

ggplot(plot.df, aes(x = label, y = rate)) +
    scale_x_discrete("Case importance [0-4]") + 
    scale_y_continuous("Proportion of cases granted permission",
                       labels = scales::percent) + 
    geom_bar(stat = "identity", colour = "grey50") +
    theme_uksc()




## ----individjudges-------------------------------------------------------
aux <- pta[,c("CASEID", "Casename", "Date", "Justice1", "Justice2", "Justice3", "Justice4", "Justice5")]
aux <- melt(aux, id.vars = c("CASEID","Casename","Date"))
aux <- subset(aux, !is.na(value))
aux$variable <- aux$value
aux$value <- 1
stopifnot(nrow(aux) == nrow(unique(aux)))

presence <- dcast(aux, CASEID + Casename + Date ~ variable, 
	value.var = "value",
	fun.aggregate = sum)

names(presence)[4:ncol(presence)] <- paste0("J",
	names(presence)[4:ncol(presence)])

### now get panel workload
wl <- merge(aux, individualWorkload,
	by.x = c("Date","variable"),
	by.y = c("Date","Judge"),
	all.x = T,
	all.y = F)

wl <- wl %>% 
	group_by(CASEID) %>% 
    summarize(panelWorkload = median(workloadTotal),
              panelWorkload.alt = min(workloadTotal))

pta <- merge(pta, wl, all.x = T, all.y = F, sort = FALSE)
pta <- merge(pta, presence, all.x = T, all.y = F, sort = FALSE)

### Now get specialisms
spec2 <- merge(aux, spec,
	by.x = "variable",
	by.y = "Judge")

spec2 <- spec2 %>% 
	group_by(CASEID) %>% 
	summarize(PanelPublic = median(Public),
		PanelFamily = median(Family),
		PanelCriminal = median(Criminal),
		PanelTax.and.Chancery = median(Tax.and.Chancery),
		PanelCivil = median(Civil),
		PanelScots = median(Scots),
		PanelNI = median(NI),
                PanelPublic.alt = max(Public),
		PanelFamily.alt = max(Family),
		PanelCriminal.alt = max(Criminal),
		PanelTax.and.Chancery.alt = max(Tax.and.Chancery),
		PanelCivil.alt = max(Civil),
		PanelScots.alt = max(Scots),
		PanelNI.alt = max(NI))

pta <- merge(pta, spec2, all.x = T, all.y = F, sort = FALSE)



## ----munge---------------------------------------------------------------

### Appellant, respondent type
pta$CentralGovt <- 0

pta$CentralGovt[which(pta$AppType %in% c("E","F","H") & !(pta$RespType %in% c("E","F","H")))] <- 1
pta$CentralGovt[which(!(pta$AppType %in% c("E","F","H")) & (pta$RespType %in% c("E","F","H")))] <- -1
                
pta$LocalGovt <- 0

pta$LocalGovt[which(pta$AppType %in% c("D","G","J") & !(pta$RespType %in% c("D","G","J")))] <- 1
pta$LocalGovt[which(!(pta$AppType %in% c("D","G","J")) & (pta$RespType %in% c("D","G","J")))] <- -1

pta$Company <- 0

pta$Company[which(pta$AppType %in% c("C") & !(pta$RespType %in% c("C")))] <- 1
pta$Company[which(!(pta$AppType %in% c("C")) & (pta$RespType %in% c("C")))] <- -1

pta$Association <- 0

pta$Association[which(pta$AppType %in% c("B") & !(pta$RespType %in% c("B")))] <- 1
pta$Association[which(!(pta$AppType %in% c("B")) & (pta$RespType %in% c("B")))] <- -1

pta$Individual <- 0

pta$Individual[which(pta$AppType %in% c("A") & !(pta$RespType %in% c("A")))] <- 1
pta$Individual[which(!(pta$AppType %in% c("A")) & (pta$RespType %in% c("A")))] <- -1

apptype.vars <- c("Individual", "Association", "Company", "LocalGovt", "CentralGovt")
test <- rowSums(pta[,apptype.vars])

if (any(test != 0)) {
    stop("Unmatched appellant/respondent pairing")
}

pta$Leapfrog <- as.numeric(grepl("EWHC|NIQB", pta$Neutral))

### First, add on one to the number of judges who heard the case
### if we have information on the first-instance court
### We add one because this is the modal number in the High Court
pta$nJudges.appord <- ifelse(pta$Leapfrog,
                             pta$nJudges,
                             pta$nJudges + 1)

### Second, add on one to the number of judges who "dissented"
### if we know that the appeal court overturned the initial decision
overturned <- grep("A", pta$APPORDDISP)
pta$nJudges.dissent[overturned] <- pta$nJudges.dissent[overturned] + 1

pta$OpinionBelow <- (1 + pta$nJudges.dissent) / (2 + pta$nJudges)



## ----areaoflaw-----------------------------------------------------------
pta$Area[which(pta$Area == "")] <- NA
pta$Area <- factor(as.character(pta$Area))
pta$Area <- relevel(pta$Area, "Civil")

pta$PanelByArea <- NA
pta$PanelByArea[which(pta$Area == "Public")] <- pta$PanelPublic[which(pta$Area == "Public")]
pta$PanelByArea[which(pta$Area == "Civil")] <- pta$PanelCivil[which(pta$Area == "Civil")]
pta$PanelByArea[which(pta$Area == "Family")] <- pta$PanelFamily[which(pta$Area == "Family")]
pta$PanelByArea[which(pta$Area == "Criminal")] <- pta$PanelCriminal[which(pta$Area == "Criminal")]
pta$PanelByArea[which(pta$Area == "Tax and Chancery")] <- pta$PanelTax.and.Chancery[which(pta$Area == "Tax and Chancery")]
pta$PanelByArea[which(pta$Area == "Scots")] <- pta$PanelScots[which(pta$Area == "Scots")]
pta$PanelByArea[which(pta$Area == "N Ireland")] <- pta$PanelNI[which(pta$Area == "N Ireland")]

pta$PanelByArea.alt <- NA
pta$PanelByArea.alt[which(pta$Area == "Public")] <- pta$PanelPublic.alt[which(pta$Area == "Public")]
pta$PanelByArea.alt[which(pta$Area == "Civil")] <- pta$PanelCivil.alt[which(pta$Area == "Civil")]
pta$PanelByArea.alt[which(pta$Area == "Family")] <- pta$PanelFamily.alt[which(pta$Area == "Family")]
pta$PanelByArea.alt[which(pta$Area == "Criminal")] <- pta$PanelCriminal.alt[which(pta$Area == "Criminal")]
pta$PanelByArea.alt[which(pta$Area == "Tax and Chancery")] <- pta$PanelTax.and.Chancery.alt[which(pta$Area == "Tax and Chancery")]
pta$PanelByArea.alt[which(pta$Area == "Scots")] <- pta$PanelScots.alt[which(pta$Area == "Scots")]
pta$PanelByArea.alt[which(pta$Area == "N Ireland")] <- pta$PanelNI.alt[which(pta$Area == "N Ireland")]



## ----pastcases-----------------------------------------------------------
windowLength <- 20

pta <- pta[order(pta$Date),]
pastvec <- pta$Area == "Public"
pastvec[is.na(pastvec)] <- FALSE

pta$pastPublic <- rollmean(pastvec, 
                           windowLength, align = "right", 	fill = 1/7)
pastvec <- pta$Area == "Civil"
pastvec[is.na(pastvec)] <- FALSE
pta$pastCivil <- rollmean(pastvec, 
                          windowLength, align = "right", 	fill = 1/7)
pastvec <- pta$Area == "Family"
pastvec[is.na(pastvec)] <- FALSE
pta$pastFamily <- rollmean(pastvec, 
                           windowLength, align = "right", 	fill = 1/7)
pastvec <- pta$Area == "Criminal"
pastvec[is.na(pastvec)] <- FALSE
pta$pastCriminal <- rollmean(pastvec, 
                             windowLength, align = "right", 	fill = 1/7)
pastvec <- pta$Area == "Tax and Chancery"
pastvec[is.na(pastvec)] <- FALSE
pta$pastTax.and.Chancery <- rollmean(pastvec, 
                                     windowLength, align = "right", 	fill = 1/7)
pastvec <- pta$Area == "Scots"
pastvec[is.na(pastvec)] <- FALSE
pta$pastScots <- rollmean(pastvec, 
                          windowLength, align = "right", 	fill = 1/7)
pastvec <- pta$Area == "N Ireland"
pastvec[is.na(pastvec)] <- FALSE
pta$pastNI <- rollmean(pastvec, 
                       windowLength, align = "right", 	fill = 1/7)

pta$PastArea <- NA
pta$PastArea[which(pta$Area == "Public")] <- pta$pastPublic[which(pta$Area == "Public")]
pta$PastArea[which(pta$Area == "Civil")] <- pta$pastCivil[which(pta$Area == "Civil")]
pta$PastArea[which(pta$Area == "Family")] <- pta$pastFamily[which(pta$Area == "Family")]
pta$PastArea[which(pta$Area == "Criminal")] <- pta$pastCriminal[which(pta$Area == "Criminal")]
pta$PastArea[which(pta$Area == "Tax and Chancery")] <- pta$pastTax.and.Chancery[which(pta$Area == "Tax and Chancery")]
pta$PastArea[which(pta$Area == "Scots")] <- pta$pastScots[which(pta$Area == "Scots")]
pta$PastArea[which(pta$Area == "N Ireland")] <- pta$pastNI[which(pta$Area == "N Ireland")]



## ----dealwithmissingdata, results = "hide"-------------------------------
keep.vars <- c("success", "Importance",
               "OpinionBelow",
               "CentralGovt", "LocalGovt",
               "Company", "Association",
               "Individual", "workloadTotal",
               "panelWorkload", "panelWorkload.alt",
               "PanelByArea", "PanelByArea.alt",
               "Area", "Leapfrog")
keeplist <- complete.cases(pta[,keep.vars])
## pta <- subset(pta, keeplist)

na2mean <- function(x) replace(x, is.na(x), mean(x, na.rm = TRUE))

## Two cases where we have structural replacement values

pta <- pta %>%
    mutate(OpinionBelow = coalesce(OpinionBelow, 1/2),
           Importance = coalesce(Importance, 0L))


tmp <- amelia(pta[, setdiff(keep.vars, "Company")], m = 1,
              boot.type = "none",
              noms = "Area")

pta.bak <- pta
pta <- tmp$imputations[[1]]


## ----smrystats-----------------------------------------------------------
impt.mean <- round(mean(pta$Importance, na.rm = TRUE), 1)
impt.sd <- round(sd(pta$Importance, na.rm = TRUE), 1)
impt.iqr <- round(quantile(pta$Importance, na.rm = TRUE,
                           probs = c(0.25, 0.75)), 1)

opbelow.mean <- round(mean(pta$OpinionBelow, na.rm = TRUE), 1)
opbelow.sd <- round(sd(pta$OpinionBelow, na.rm = TRUE), 2)

panel_wl.mean <- round(mean(pta$panelWorkload, na.rm = TRUE), 1)
panel_wl.sd <- round(sd(pta$panelWorkload, na.rm = TRUE), 1)

panel_wl_alt.mean <- round(mean(pta$panelWorkload.alt, na.rm = TRUE), 1)
panel_wl_alt.sd <- round(sd(pta$panelWorkload.alt, na.rm = TRUE), 1)

ct_wl.mean <- round(mean(pta$workloadTotal, na.rm = TRUE))
ct_wl.sd <- round(sd(pta$workloadTotal, na.rm = TRUE))

spec_median.mean <- round(mean(pta$PanelByArea, na.rm = TRUE), 1)
spec_median.sd <- round(sd(pta$PanelByArea, na.rm = TRUE), 1)

spec_max.mean <- round(mean(pta$PanelByArea.alt, na.rm = TRUE), 1)
spec_max.sd <- round(sd(pta$PanelByArea.alt, na.rm = TRUE), 1)



## ----mockfig, fig.cap = "Illustrative model of logistic regression coefficients. See text for a fuller description. ", out.width = "100%"----
plot.df <- data.frame(variable = LETTERS[1:5],
                      y = c(1, 0.5, -0.5,
                            0.4, -0.3),
                      se = c(0.24, .33, 0.1,
                             0.24, 0.24),
                      annot = c("Very good evidence success goes up with A",
                                "Good evidence success goes up with B",
                                "Very good evidence success goes down with C",
                                "Very good evidence that D...",
                                "... is larger than E"))

plot.df$annot <- str_wrap(plot.df$annot, 20)

plot.df$hi <- plot.df$y + qnorm(39/40) * plot.df$se
plot.df$lo <- plot.df$y - qnorm(39/40) * plot.df$se

plot.df$hhi <- plot.df$y + qnorm(11/12) * plot.df$se
plot.df$llo <- plot.df$y - qnorm(11/12) * plot.df$se

plot.df$isSignificant <- apply(plot.df[,c("lo", "hi")], 1, isSig)

inner.bar.col <- grey(0.7)
outer.bar.col <- grey(0.5)

## Red-y ones
inner.bar.col <- "#fcae91"
outer.bar.col <- "#a50f15"

inner.bar.size <- 0.5
outer.bar.size <- 1

plot.df$inner.bar.col <- ifelse(plot.df$isSignificant,
                                inner.bar.col,
                                "lightgrey")
plot.df$outer.bar.col <- ifelse(plot.df$isSignificant,
                                outer.bar.col,
                                "darkgrey")

my.ylims <- range(c(plot.df$lo, plot.df$hi))
## Over-ride
my.ylims[2] <- 1.8

sec.ylab <- "Odds are ... times greater"

plot.df$variable <- factor(plot.df$variable,
                           levels = rev(LETTERS[1:5]),
                           ordered = TRUE)

retval <- ggplot(plot.df,
                 aes(x = variable,
                     y = y,
                     ymin = lo,
                     ymax = hi)) +
    scale_x_discrete("") +
    scale_y_continuous("Coefficient value",
                       limits = my.ylims,
                       sec.axis = sec_axis(~exp(.),
                                           name = sec.ylab)) +
    geom_linerange(aes(colour = inner.bar.col),
                       size=inner.bar.size, 
                   alpha = 0.8) +
    geom_linerange(aes(ymin = llo,
                       ymax = hhi,
                       colour = outer.bar.col),
                   size=outer.bar.size, 
                   alpha = 0.8) +
    geom_hline(aes(yintercept=0), lty=1, colour = "darkgrey") +
    geom_point(aes(colour = outer.bar.col),
               size = (inner.bar.size + outer.bar.size) * 2.5,
               shape = my.pchs[1],
                   stroke = 1.5) +
    geom_text(aes(x = variable, y = hi,
label = annot, hjust = 0), size = 1.5) + 
    scale_colour_identity() +
    scale_fill_identity() + 
    theme_uksc() +
    coord_flip()

print(retval)



## ----legalmodfig, fig.cap = "Legal model of permission to appeal decisions. Thick lines indicate 83 percent confidence intervals. Thin lines indicate 95 percent confidence intervals. ", fig.scap = "Legal model of permission to appeal decisions", fig.width=7, fig.height=4----

## Figure height should be what, (nCoefs + 2) / 2.5?
pta$Area <- as.character(pta$Area)
pta$Area[is.na(pta$Area)] <- "Unclear"
pta$Area <- factor(pta$Area)
pta$Area <- relevel(pta$Area, "Civil")

pta$Importance <- as.double(pta$Importance)

t_prior <- student_t(df = 7, location = 0, scale = 2.5)

if (file.exists("cache/03_cached_results.RData")) {
    load("cache/03_cached_results.RData")
    } else { 
        legalmod <- bayesglm(success ~ rescale(Importance) +
                                 rescale(OpinionBelow) +
                                 Area,
                             data = pta,
                             family = binomial)


        legalmod.rstanarm <- stan_glm(success ~ rescale(Importance) + rescale(OpinionBelow) + Area,
                                      data = pta,
                                      family = binomial(link = "logit"), 
                                      prior = t_prior, prior_intercept = t_prior,  
                                      chains = 2, cores = 2, seed = 1982, iter = 1000,
                                      refresh = -1, verbose = FALSE)

        org_mod <- bayesglm(success ~ rescale(panelWorkload) + rescale(workloadTotal) + rescale(PanelByArea),
                            data = pta,
                            family = binomial)

        org_mod.alt <- bayesglm(success ~ rescale(panelWorkload.alt) + rescale(workloadTotal) + rescale(PanelByArea.alt),
                                data = pta,
                                family = binomial)

        org_mod.rstanarm <- stan_glm(success ~ rescale(panelWorkload) +
                                         rescale(workloadTotal) + rescale(PanelByArea),
                                     data = pta,
                                     family = binomial(link = "logit"), 
                                     prior = t_prior, prior_intercept = t_prior,  
                                     chains = 2, cores = 2, seed = 1982, iter = 1000,
                                     refresh = -1, verbose = FALSE)

        org_mod.alt.rstanarm <- stan_glm(success ~ rescale(panelWorkload.alt) +
                                             rescale(workloadTotal) + rescale(PanelByArea.alt),
                                         data = pta,
                                         family = binomial(link = "logit"), 
                                         prior = t_prior, prior_intercept = t_prior,  
                                         chains = 2, cores = 2, seed = 1982, iter = 1000,
                                         refresh = -1, verbose = FALSE)


        polmod <- glm(success ~ CentralGovt + LocalGovt +
                                        # Company +
                          Association + Individual +
                          Area,
                      data = pta,
                      family = binomial)

        polmod.rstanarm <- stan_glm(success ~ CentralGovt + LocalGovt +
                                        # Company +
                                        Association + Individual +
                                        Area,
                                    data = pta,
                                    family = binomial(link = "logit"), 
                                    prior = t_prior, prior_intercept = t_prior,  
                                    chains = 2, cores = 2, seed = 1982, iter = 1000,
                                    refresh = -1, verbose = FALSE)

        combmod <- bayesglm(success ~ rescale(Importance) + rescale(OpinionBelow) + Area +
                                rescale(panelWorkload) + rescale(workloadTotal) + rescale(PanelByArea) +
                                CentralGovt + LocalGovt +
                                Association + Individual,
                            data = pta,
                            family = binomial)

        combmod.alt <- bayesglm(success ~ rescale(Importance) + rescale(OpinionBelow) + Area +
                                    rescale(panelWorkload.alt) + rescale(workloadTotal) + rescale(PanelByArea.alt) +
                                    CentralGovt + LocalGovt +
                                    Association + Individual,
                                data = pta,
                                family = binomial)

        combmod.rstanarm <- stan_glm(success ~ rescale(Importance) + rescale(OpinionBelow) + Area +
                                         rescale(panelWorkload) + rescale(workloadTotal) + rescale(PanelByArea) +
                                         CentralGovt + LocalGovt +
                                         Association + Individual,
                                     data = pta,
                                     family = binomial(link = "logit"), 
                                     prior = t_prior, prior_intercept = t_prior,  
                                     chains = 2, cores = 2, seed = 1982, iter = 1000,
                                     refresh = -1, verbose = FALSE)
        save(legalmod.rstanarm, org_mod.rstanarm, org_mod.alt.rstanarm,
             polmod.rstanarm, combmod.rstanarm,
             file = "cache/03_cached_results.RData")
    }

print(mycoefplot(legalmod.rstanarm,
           ylab = "Coefficient value",
           sec.ylab = "Odds are ... times greater",
           variable.labels = list("(Intercept)" = "Intercept",
                                  "rescale(Importance)" = "Importance",
                                  "rescale(OpinionBelow)" = "Opinion below",
                                  "AreaCriminal"='Area: Criminal',
                                  "AreaFamily"= 'Area: Family',
                                  "AreaN Ireland"= 'Area: N Ireland',
                                  "AreaPublic" = 'Area: Public',
                                  "AreaScots" = 'Area: Scots',
                                  "AreaTax and Chancery" ='Area: Tax and Chancery')))
        


## ----pulllegal-----------------------------------------------------------
impt.coef <- tidy(legalmod.rstanarm) %>%
    filter(term == "rescale(Importance)") %>%
    dplyr::select(estimate) %>%
    pull() %>%
    round(digits = 2)
impt.exp <- round(exp(impt.coef), 2)

scots.coef <- tidy(legalmod.rstanarm) %>%
    filter(term == "AreaScots") %>%
    dplyr::select(estimate) %>%
    pull() %>%
    round(digits = 2)
scots.exp <- round(exp(scots.coef), 2)



## ----orgmods, fig = TRUE, fig.cap = "Organisational model of permission to appeal decisions. Thick lines indicate 83 percent confidence intervals. Thin lines indicate 95 percent confidence intervals. Coefficients have been standardised. ", fig.width = 6, fig.height = 6----
## Figure height should ordinarily be (nCoefs + 2) / 2.5
## But we have two models, so it's more than that 
### Add in custom ylims
## my.ylims <- c(NA, NA)
## my.ylims[1] <- min(c(tidy(org_mod.rstanarm)[,2] - 2 * tidy(org_mod.rstanarm)[,3],
##       tidy(org_mod.alt)[,2] - 2 * tidy(org_mod.alt)[,3]))
## my.ylims[2] <- max(c(tidy(org_mod)[,2] + 2 * tidy(org_mod)[,3],
##       tidy(org_mod.alt)[,2] + 2 * tidy(org_mod.alt)[,3]))

p1 <- mycoefplot.mult(modlist = list(org_mod.rstanarm, org_mod.alt.rstanarm),
                      model.names = c("Majority rule", "Rule-of-one"),
                      ## ylim = my.ylims,
                      ylab = "Coefficient value",
                      sec.ylab = "Odds are ... times greater",
                      variable.labels = list("(Intercept)" = 'Intercept',
                                        "rescale(panelWorkload)" = 'Panel workload',
                                        "rescale(workloadTotal)" = 'Court workload',
                                        "rescale(PanelByArea)" = 'Panel specialism in area',
                                        "rescale(panelWorkload.alt)" = 'Panel workload',
                                        "rescale(PanelByArea.alt)" = 'Panel specialism in area'))

print(p1)



## ----pullpolitical-------------------------------------------------------
cgov.coef <- tidy(polmod.rstanarm) %>%
    filter(term == "CentralGovt") %>%
    dplyr::select(estimate) %>%
    pull() %>%
    round(digits = 2)
cgov.exp <- round(exp(cgov.coef), 2)



## ----polmodout, fig = TRUE, fig.cap = "Political models of permission to appeal decisions. Thick lines indicate 83 percent confidence intervals. Thin lines indicate 95 percent confidence intervals. Coefficients have been standardised. ", fig.width = 6, fig.height = 4.8----
mycoefplot(polmod.rstanarm,
           ylab = "Coefficient value",
           sec.ylab = "Odds are ... times greater",
           variable.labels = list(`CentralGovt` = 'Appellant is central govt.',
                                  `LocalGovt` = 'Appellant is local govt.',
                                  `Association` = 'Appellant is an association',
                                  `Individual` = 'Appellant is an individual',
                                  `AreaCriminal` = 'Area: Criminal',
                                  `AreaFamily` = 'Area: Family',
                                  `AreaN Ireland` = 'Area: N Ireland',
                                  `AreaPublic` = 'Area: Public law',
                                  `AreaScots` = 'Area: Scots',
                                  `AreaTax and Chancery` = 'Area: Tax and Chancery'))



## ----combmod,fig.width=6, fig.height=6.8, fig = TRUE, fig.cap = "Combined model of permission to appeal decisions, majority rule. Thick lines indicate 83 percent confidence intervals. Thin lines indicate 95 percent confidence intervals. ", out.width = "100%"----
## Fig height for one would be = (15 + 2) / 2.5 = 6.8

p1 <- mycoefplot(combmod.rstanarm,
                 ylab = "Coefficient value",
                 sec.ylab = "Odds are ... times greater",
                 variable.labels = list("rescale(Importance)" = 'Importance',
                                        "rescale(OpinionBelow)" = 'Opinion below',
                                        "rescale(panelWorkload)" = 'Workload, panel',
                                        "rescale(panelWorkload.alt)" = 'Workload, panel',
                                        "rescale(workloadTotal)" = 'Workload, court',
                                        "rescale(PanelByArea)" = 'Specialisation',
                                        "rescale(PanelByArea.alt)" = 'Specialisation',
                                        "CentralGovt" = 'Appellant is central govt.',
                                        "LocalGovt" = 'Appellant is local govt.',
                                        "Association" = 'Appellant is an association',
                                        "Individual" = 'Appellant is an individual',
                                        "AreaCriminal" = 'Area: Criminal',
                                        "AreaFamily" = 'Area: Family',
                                        "AreaN Ireland" = 'Area: N Ireland',
                                        "AreaPublic" = 'Area: Public law',
                                        "AreaScots" = 'Area: Scots',
                                        "AreaTax and Chancery" = 'Area: Tax and Chancery'))

if (FALSE) { 
p2 <- mycoefplot(combmod.alt,
                 ylab = "Coefficient value",
                 sec.ylab = "Odds are ... times greater",
                 variable.labels = list("rescale(Importance)" = 'Importance',
                                        "rescale(OpinionBelow)" = 'Opinion below',
                                        "rescale(panelWorkload)" = 'Workload, panel',
                                        "rescale(panelWorkload.alt)" = 'Workload, panel',
                                        "rescale(workloadTotal)" = 'Workload, court',
                                        "rescale(PanelByArea)" = 'Specialisation',
                                        "rescale(PanelByArea.alt)" = 'Specialisation',
                                        "CentralGovt" = 'Appellant is central govt.',
                                        "LocalGovt" = 'Appellant is local govt.',
                                        "Association" = 'Appellant is an association',
                                        "Individual" = 'Appellant is an individual',
                                        "AreaCriminal" = 'Area: Criminal',
                                        "AreaFamily" = 'Area: Family',
                                        "AreaN Ireland" = 'Area: N Ireland',
                                        "AreaPublic" = 'Area: Public law',
                                        "AreaScots" = 'Area: Scots',
                                        "AreaTax and Chancery" = 'Area: Tax and Chancery'))
}

print(p1)



## ----pullcomb------------------------------------------------------------
publaw.coef <- tidy(combmod.rstanarm) %>%
    filter(term == "AreaPublic") %>%
    dplyr::select(estimate) %>%
    pull() %>%
    round(digits = 2)
publaw.exp <- round(exp(publaw.coef), 2)



## ----fitstats, fig = TRUE, fig.cap = "Goodness-of-fit statistics for different models. Thin lines indicate 95 percent credible intervals; thick lines indicate 83 percent credible intervals. ", warning = FALSE----
epcp.rstanarm <- function(mod) {
    poslp <- posterior_linpred(mod)
    preds <- (plogis(poslp) > 0.5)
    
    actual <- matrix(mod$y,
                     nrow = 1000,
                     ncol = length(mod$y),
                     byrow = TRUE)
    correct <- (preds == actual)

    correct <- apply(correct, 1, mean)
    epcp <- data.frame(term = "PCP\n(bigger is better)",
                       mean = mean(correct),
                       lo = quantile(correct, 0.025),
                       hi = quantile(correct, 0.975),
                       llo = quantile(correct, 1/12),
                       hhi = quantile(correct, 11/12))

    the.loo <- loo(mod)
    looic <- the.loo$estimates["looic", 1]
    loose <- the.loo$estimates["looic", 2]
    the.loo <- data.frame(term = "LOOIC\n(smaller is better)",
                          mean = looic,
                          lo = looic + qnorm(1/40) * loose,
                          hi = looic + qnorm(39/40) * loose,
                          llo = looic + qnorm(1/24) * loose,
                          hhi = looic + qnorm(23/24) * loose)
    return(rbind(epcp, the.loo))
}

if (file.exists("cache/03_gof.RData")) {
    load("cache/03_gof.RData")
} else { 

    legal.stats <- epcp.rstanarm(legalmod.rstanarm)
    org.stats <- epcp.rstanarm(org_mod.rstanarm)
    pol.stats <- epcp.rstanarm(polmod.rstanarm)
    comb.stats <- epcp.rstanarm(combmod.rstanarm)

    legal.stats$Model <- "Legal"
    org.stats$Model <- "Organisational"
    pol.stats$Model <- "Political"
    comb.stats$Model <- "Combined"

    null.stats <- data.frame(Model = rep("Null", 2),
                             term = c("PCP\n(bigger is better)",
                                      "LOOIC\n(smaller is better)"),
                             mean = c(1 - mean(legalmod.rstanarm$y,
                                               na.rm = TRUE),
                                      NA))

    gof.df <- rbind(
        rbind(
            rbind(legal.stats,
                  org.stats),
            pol.stats),
        comb.stats)

    gof.df <- merge(gof.df, null.stats, all = TRUE)

    my.levels <- gof.df %>%
        filter(grepl("LOOIC", term)) %>%
        arrange(mean) %>%
        pull(Model)

    gof.df$Model <- factor(gof.df$Model,
                           levels = rev(my.levels),
                           ordered = TRUE)

legal.loo <- loo(legalmod.rstanarm)
combined.loo <- loo(combmod.rstanarm)
political.loo <- loo(polmod.rstanarm)
org.loo <- loo(org_mod.rstanarm)

    save(list = c("legal.loo", "combined.loo",
                  "political.loo", "org.loo",
"gof.df"),
         file = "cache/03_gof.RData")

}

    ggplot(gof.df, aes(x = Model, y = mean, ymax = hi, ymin = lo)) +
        geom_linerange(size = inner.bar.size,
                       colour = inner.bar.col,
                   alpha = 0.8) +
        geom_linerange(aes(ymin = llo, ymax = hhi),
                       size = outer.bar.size,
                       colour = outer.bar.col,
                       alpha = 0.8) +
    geom_point(aes(colour = outer.bar.col),
               size= (inner.bar.size + outer.bar.size) * 2.5,
               shape= my.pchs[1],
               stroke = 1.5) +
    scale_colour_identity() + 
    facet_wrap(~term, scales = "free_x") +
    scale_y_continuous("Value") +
    coord_flip() +
    theme_uksc()
